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ABSTRACT 


The objective of this work was to analyze a computer 
program using three dimensional quadratic isoparametric 
finite elements for structural analysis. Three problems 
with classical solutions were run with various mesh sizes 
Nine tae computer program being tested. The data computed 
was then extensively analyzed, and compared with the 
РЕЙСЕ БСА solutions. The analysis of a fourth problem was 
continued and compared with results obtained in an earlier 


project. 
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1. INTRODUCTION 


Before the advent of the finite element technique and 
рпе digital computer, the solutions T non-trivial 
elasto-static problems were unobtainable. Since that time 
many finite elements have been devised with corresponding 
si er programs for their use. 

One of the most versatile family of finite elements was 
developed by Professor Zienkiewicz and his co-workers at 
the University of Wales, Swansea, U.K. This family of 
elements, called "ISOPARAMETRIC", has been investigated at 
une Naval Postgraduate School by Professor G. Cantin, and a 
stress analysis program called TRISOP was written using a 
¿0 nodal point quadratic element. 

Four problems with known solutions will be solved using 
ШЕЕСОР апа analyzed in this thesis. From this analysis 
conclusions will be made as to the effectiveness and accuracy 


PETRISOR: 





ИШ МОВІ ЗОРАНАМЕТЕВНІС ЕТМІТЕ ЕБЕМЕМТ5 


i S l C Co Tonal element discussed in this chapter 
wache type of element used in the program analyzed by this 
MD Mg = пиле made to fully describe this 
ype of element, but just to give some insight as to how 
ane element is constructed and designed. 

The three-dimensional quadratic isoparameteric finite 
element, hereafter referred to as the "element", is a cube 
ПЕ ат! Sides two units in length in іс undeformed, non- 
dimensional state. The element has 20 nodal points and is 
oriented with its geometric center at the origin of its 
non-dimensionalized coordinates, (&,n,¢), and transferred 


to global, (X,Y,2), coordinates as shown in Figure 1. 


A. TRANSFORMATIONS 
tranc кеа from global coordinates to non- 
dimensional coordinates are accomplished using shape functions 


ον ΘΠΟΠΠ τα equations 1 


О = N, (E,n,z)x, 
СЕ τ. (1) 
2(Е,п,6) = N, (&,n,z)z, 


where Ns (£,n, t) are shape functions, and x,,y,,2, are 


coordinates of the nodal points. Similarly. displacement 


transformations are made using equations 2, 








Figure 1. Element Transformation 


10, 





πμ... N. (&,n,z)u, 


cem Ee) " N, (E,n,z)v; (2) 


ν(ξ of 2c) N, (š ЭТ 5S ји. 


where u, v, and w are displacements in the global (X,Y,Z) 
ше erence frame. 

DN s Or C lon or a volumetric increment from 
loeal (ξ,η.ζ) to the global (X,Y,Z) system of reference is 


shown in equation 3. The Jacobian, [J] used in the 
dxdydz = det[J] d&dndz (22?) 


transformation is shown in equation 4. 


ax ay az 
DIG 9 š οξ 
- | 3x әу az 
[J] -: - ~ (4) 
ax ay az 
οζ 9 Ç 95 


Б: ӘНДЕР FUNCTIONS 

For an element with nodes numbered as in Figure 2, the 
shape functions obtained from Reference 1 are shown in 
Сема тоне. 5. 


ТООС осе 1: 32. 5, 7T. 13, 15, 17, апа 19 


ЕК eee τει επι Ἑ -2) (5) 


where ES = TE and Es се ноте огју for Me and ζω 
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Figure 2. Element Nodal Point Numbering System 


T2 





Midside Nodes: 2, 6, 14, and 18 
- 2 
ы (1/4)(1 - €°)(1 + Е + 55) 


Midside Nodes: 4, 8, 16, and 20 (5) 
Е 2 
ο... ασ πα περα τς) 
Midside Nodes: 9. 10, 11, апа 12 
= 2 
N, = (1/0) (1 – с) + ЕЈ + n.) 
Опе should take note that if any other nodal numbering 
Er Zen wich respect со E,n,c is used, the nodal point 
numbers that go with the above shape functions must be 


аси сеа accordingly. It should also be noted that the 


program under analysis uses the system mentioned above. 


C. COMPUTING STRESS AND STRAIN 


Noain? strains are computed using the following 


relationship. 
9 
ες το Ü | 
д 
кү 0 зу 0 
u 
9 
Ez 0 E T 
- V (6) 
9 9 
Үху ду 9X : η 
9 9 
Yxz 97 0 9 X 
9 9 
Yyz Б 22 3y 


where u, v, and w are obtained from equation 2. The deri- 


vatives of the displacements in global X, Y, Z coordinates 


απο piven by: 
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а? 
C 
а 
« 
Ω7 
>, 
а? 
c 
а? 
< 
а? 
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9 Χ ðX 9 X οξ οξ οξ 
ди о“ ом Т -1 ди 9 V дм 
ду ду ду LJ] 9η an 3n (7) 
au ὃν ам au ὃν Әм 
97 д2 92 95 9 G ος 


where [J] is the Jacobian matrix (equation 4). The stresses 


are computed from strains by: 


- At2G PA À 0 0 0 3% 
σν λ αλ 0 0 0 Εν 
ο λ А A+2G 0 0 0 e, | (8) 
"Ху I 0 0 0 G 0 0 Lum 
Xo 0 0 00 0 G 0 ЕШ 
Туя 0 0 0 0 0 G Yyz 
= = i (9) 


in which: G = any > À * (yoo 


5 modulus, and vw is Poisson's ratio. 
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ISO EA QUADRATTC ISOPARAMETRIC 
FINITE ELEMENT PROGRAM 
Πο ο лані  реююртшуе of this work is to evaluate TRISOP, 
and give some insight into its practical use. Reference 1 
КОШ ы ыле theory and methods used in the program, and 


тиште 3 15 а simplified flow diagram. 


A. INPUT DATA 
The input data required by TRISOP are: 


The number of elements 


The total number of nodal points 


The number of different materials 


The block size for the large capacity solver 
Шет Оет of nodal peints with boundary conditions 
The number of nodal points with concentrated loads 


Element connectivity for each element 


The coordinates of each nodal point 


ο CON NUT + GU) n5 HF 


Young's Modulus, Poisson's Ratio for each material 


= 
о 


Concentrated loads 


mee Boundary conditions 
The element connectivity is a correlation between the over- 
all mesh numbering system and the standard numbering system 
ecc tes nte пете. 

When large problems are to be solved, preparing the 
connectivity and coordinate input is not a trivial task. 
Also, if there is a pressure or gravity load, the computa- 
tions needed to convert into corresponding nodal point 
toads can pe extremely time consuming. There is, however, 


a mesh generator program [2] that will solve this problem. 


1» 
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ESE 


| MERGE ο 


FLOAD (2) 





= — > = 


кш πο | (2) 
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Figure 3. Functional Flow Diagram 
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With a minimum input, the mesh generator will compute the 
υπ ιν eonmeecivicy, nodal point coordinates with input in 
DeAnna ical or rectangular coordinates, and nodal point 
loads for pressure or gravity loads. The mesh generator 
will print the output, draw a two dimensional picture of 
ШІС πε5η and punch data cards for TRISOP. The only other 
input to TRISOP that can require more than a few cards is 
mes boundary conditions. The mesh generator will not 
Ν.Ε DOUNnGary condition cards, but they are easily and 


Q hly produced by hand. 


В. TECHNIQUES FOR MOST EFFICIENT USE OF TRISOP 

ПИ КО пио S presenti orm, has a constant core storage 
requirement. The two variables with the size of the problem 
КШ 1555 5оогаре requirements and running time. To make the 
most efficient use of TRISOP, it is essential to take advan- 
Cage of symmetry whenever possible. It is also essential 
to design the problem mesh in order to reduce the half-band 
width c? the resulting system of equations to a minimum. 

ЕЕЕ Е ааа width is a function of the difference 
between the highest and lowest nodal point number in any 
element. To obtain the smallest band width one must start 
numbering on the face having the least number of elements 
along the side having the smallest number of elements as 
shown in Figure Y. 

TRISOP does not give exact answers to a problem, but 


will converge asymptotically and monotonically to an exact 


nq 
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Figure 4. 


Numbering System for Minimum Half-Band Width 
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БЕГІС ИО ГЕШ ВОМ mesh refinement. To insure reliable 
results to a problem, a convergence study must be made. 
The convergence technique used here consists of plotting 
the displacement of a convenient nodal point, versus 1/N* 
frre N is the number of elements in the mesh. If three 
та prot in a straight line, the extrapolation to the 
Ex is Justified. This convergence technique is an 
adaptation of a technique developed by L. P. Richardson [3]. 
Поло лес о о 5 technique Was developed for extrapolating 
the results of central finite difference approximations 
where the truncation error is of the order (Ж and h is 
Цитен осе difference interval. When the truncation error 


is of the order (h°) the extrapolated value is found by: 


s Ши. 
A xtrap B m 2 (10) 


lum formula as plotted in Figure 5. 





Figure 5. Richardson's Technique 
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шаға third value, x oma with a third interval, h 


5 3° 
and х3 plots on the extrapolation line, the extrapolated 
value is presumed to be exact [3]. 

There is no guarantee that the error in TRISOP is of 
order 1/N*, but this method of ascertaining convergence 
was adopted after some numerical experimentation. Since 
it is assumed that the solution approaches the exact solu- 
tion asymptotically, if three deflections plot in a straight 
ШЕСЕ 15 assumed that extrapolation is valid. If they 
don't plot in a straight line there is no way of predicting 
the accuracy of the extrapolated result. This technique 


has worked quite well with many problems, but a better 


means of determining convergence is needed. 


DAR TMPROVEMENTS MADE TO TRISOP 

Gaussian integration is used extensively in the compu- 
Шол Of the element stiffness. Initially, four Gauss 
pounes in each of the three directions of a coordinate 
system were used in the solution. The subroutine in which 
this was accomplished was called CUBA. During a visit by 
Ἱ ов" о C. Zienkiewicz this author was told that using 
two Gauss points in the integration improved the solution. 
ΠΕ ΒΘ, а two point integration subroutine, was 
КОШЕ КОШЕ tor CUBA in TRISOP. This change yielded better 
results with a coarser mesh for all problems tested. The 
change from CUBY to CUB2 also reduces significantly the 
το CPU time for the calculation of element 


stiffness. 
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If an element is deformed into an extreme shape such 
as the element in Figure 6, the Jacobian becomes singular 
ο οπής similar to 3, 10, and 15. Since equation 7 
which uses the inverse of the Jacobian is used in computing 
stress and strain values, the computer algorithm fails at 
ШіШ:е 2, 10, and 15. To eliminate this singularity, 
TRISOP was modified to displace each node by a small dis- 
Lance away from its actual location at the time the 
Jacobian is formed for that node. The flow of computations 
is then uninterrupted, and the results for all other nodes 
are unaffected. The stresses and strains computed for nodes 
witha singularity are meaningless and should be eliminated 


lNusmecny further considerations. 


2L 








Degenerate Element 


Figure 6. 
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IV. SOLUTIONS TO CLASSICAL PROBLEMS USING TRISOP 


Four classical problems were considered in the analysis 
of TRISOP. The four problems were a simply supported bean, 
a pinched disk, the Boussinesq problem, and a pinched 


ey kinder. 


ПР» ІМРІУ SUPPORTED BEAM 
The simply supported beam shown in Figure 7 was analyzed. 
IO ta ssical Solution 
The classical solution used was an Airy stress 
function developed using elasticity theory [4]. The Airy 
stress function for the problem under consideration is 


shown in equation 11. 


2 
й 2. 39 2 2 3 а 
фе - £ x + бе ху - уху Шыл 5) ν ην (11) 


Emustion lr satisfies equation 12 as is required 


Y Ц Ц 
Eg Ue I 0 (12) 
9X әх әу ду 


for a valid Airy stress function in the absence of body 
forces. Stress values are found from equation 11 using 


equation 13. 





E Š . = 3 $ ; ЖА - т сло 
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Evaluating these relations yields equation 1l 


Beo, atu. Ee 
sr x ( 5 y? - ey + = c2) (14) 
i, із (ES - y*)x 

where I = — is the moment of inertia of the beam. Since 


Шоссе по body forces, such as gravity loading, this is 
the solution to both the plane stress and plane strain 
problems. The stress component, sigma xis an exact solu- 
tion to the problem only if the distributed normal force 


shown in equation 15 is applied 


:-. ne Ў E ( 5 y? = 


σ 


мл го 


б) | CEDE 


Пасо епалог сре beam. Since this force distribution 
Dass zero resultant force, and zero resultant БЕА. ic can 
be concluded by Saint-Venant's Principle [4] that sigma x 
ШОО ш оше distance from the ends of the beam. 
saint -Venant's Principle assumes that localized forces іп 
static equilibrium (it лен ос sro localized stresses 
audsstrains. 

The displacement of the center of the beam is 


found using equation 16. This displacement is greater 


no 


- MC Ц ν 
ο. S [1 + RES 2 (= + 5) ] (16) 





than the displacement found using elementary theory because 
mere the assumption is made that cross sections of the beam 
ecmain plane during bending. 
2. TRISOP Formulation 

Considering the problem in Figure 7, there is a 
plane of symmetry parallel to the XY plane passing through 
Meno axis at .5 inches. All the points in this plane of 
Eee cry were given a geometrical constraint which set the 
dısplacement component in the Z direction equal to zero. 
ο... μαι ОГ the beam from 0 < Z < „5 was considered. For 
mmemmodal points in the XZ plane at X = tL, the displace- 
ment component in the Y direction was set equal to zero. 
Pally, Tor the nodal points on the YZ plane the displace- 
l component in the X direction was set equal to zero. 
The distributed load of 1000 psi is modeled using a consist- 
ООШ Сас vector [5]. This vector dictates that the force 
on the face of each loaded element be divided as shown in 


Figure 8. 





Figure 8. Consistent Load Vector 





Шо Results 

A cOnvergence study of the deflections at the 
ЕСШЕ Of the beam is shown in Figure 9 using the meshes 
shown in Figure 10. The results indicate that extrapola- 
tion to a mesh composed of an infinity of elements is 
Justified. The classical value for the deflection at the 
center of the beam is 0.0180 inches. The extrapolated 
value given by TRISOP was 0.0188 inches, avalue .0008 
memes greater than the classical value. The classical 
solution, being a two dimensional solution, assumes that 
штете is no deflection perpendicular to the XY plane. 
Шо Makes the classical solution stiffer than the three 
dimensional solution. Figure 11 shows the displacement 
Ewnponent іп the Z direction in the XY plane given by TRISOF. 

When TRISOP results for sigma x, sigma y, and tau xy 
Lax conpared with the classical solution using nodes on 
Miemmesm сасе the values given by sigma x were within at 
least .1 percent of the classical value τος all nodal 
Ponts. Оп the other hand, the values for sigma y and 
tau xy were not nearly as good. It was noted, however, Lia 
the best results were obtained from interior nodes in the 
mesh. А 8х3х2 mesh was then analyzed using data from the 
mid plane parallel to the XY plane. The results from this 
anal ss are shown in Figures 12, 13, and 14. Since the 
beam is symmetrical about the origin along the X axis, 


suu eft end of the beam, -L « X « 0, is shown. 
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8x2x1 mesh 





8x3x1 mesh 





8x4x1 mesh 


Figure 10. Simply Supported Beam Meshes 
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Displacements x 1072 on XY Plane of the 
Simplyeoupperved Beam- 


Figure ll 
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Guasscical Solution 
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X=-L Classical Solution X=0 
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Cfassical Solution Х-0 


24 20 16 12 8 i 0 


The lower half is the image of the upper 
half. The right end is the image of the 
left end with reversed signs. 


ү πι ος о поз оп X=0 


24 20 16 12 8 4 


Figure 14. Simply Supported Beam 
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ц. Conclusions 
The results of this analysis indicate that much 
more reliable data is obtained from. interior nodes. The 
poor results obtained for sigma y and tau xy on the mesh faces 
us Ἱπαισοαστε that the values obtained under a consistent 
load vector are better in planes where the loaded nodal 
points have four common elements as in the mid plane of the 


8x3x2 mesh. 


ЕЙ PINCHED DISK 

A pinched disk with two equal and opposite forces 
acting on the diameter as shown in Figure 15 was analyzed. 

lec lLassical solution 

А сйаввтсаі two dimensional solution was devised 

using an Airy stress function by H. Hertz, and is discussed 
in detail in Reference 4. The Airy stress function, 
equation 17, gives the stress in terms of 6 and r for each 


load as shown in Figure 16. The radial stress, the only 
P 
Бо = r6 sin 9 NT) 
non zero stress, is shown in Equation 18. Equation 18 


2 
кро ба at) 


а 
її 
КЕ” 
о? 
= 
4 


Q2 
H 
3 IP 


leads to the presence of an isotropic state of compression 
Поета а су Рута all around the radial surface of the 
ОШ КЛ СЫ рос the Doundary of this unwanted stress an 


isotropic tension equal to 2Ρ/πᾶ is added to the disk. 
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Figure 16. Pinched Disk Radial Stress 


36 








WEC Converting sigma τ and sigma ro WC ony coordinates, 


and adding the stresses from each load together, equations 


19 result. 
σ = —0. sin’ 8. - c zum 6. + er 
X r 1 р 2 та 
Лі 2 
σ = y BOSS 8. = 0 soe 6. + er (19) 
y р 18 г 2 тд 
T 2 
“xy == ο θη σος θη + pn 9, cos 9, 


ΙΕ πορππ ασ σα 
The problem was formulated for TRISOP as shown іп 
Ерге wl, using half спе thickness of the disk. The 
boundary conditions were w = 0 in the XY plane, u = Din 
the YZ plane, and v = 0 in the XZ plane. Figure 18 shows 
one face of the 8x8xl mesh. To model the load, a consistent 
load vector was again used [5] with values as shown in 
Ligure 17. 
3. Results 
The deflections on the radial surface of the disk 
at a distance of 0.975 inches from the load were 0.01130 
inches for the 4x4x1 mesh, and 0.01133 inches for the 
8x8xl mesh. Due to the small difference in these two 
deflections, it was assumed that the solution could be 
Миа оне е с оно еп icant error. The convergence 


БИНЕШ слома in Figure 19. 
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ШОШ ο ὉἙνποισα Disk, TRISOP Formulation 
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Pinched Disk 8x8x1 Mesh 
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kieure 19. Pinched Disk Convergence Study 


40 





Contour graphs were drawn using data from the 8x8x1 
mesh to compare TRISOP data with the classical solution 
obtained by Carlos Felippa [6]. Two graphs were drawn, one 
Morag πια side node data, and one using corner node data 
ПЕЕЛИ Спе results shown in Figure 20. The plots for sigma y 
and tau xy using corner and mid side nodes are shown in 
heures 21 and 22. All TRISOP plots shown used data from 
the constrained (ХУ) plane. However, the data on the free 
plane (Z = .5 inches), and the mid plane (Z = .25 inches) 
was virtually the same as that in the corresponding nodal 
ЕСІП 5 in the XY plane. 

4. Conclusions 

The TRISOP solution agrees with the classical solu- 
tion for sigma y, tau xy, and the mid side node data for 
Sigma x. The results for sigma x using corner node data 
indicates a much higher state of stress in the X direction 
МЕЕ ОГУ Under the load. It should be noted, however, that 
there are only two nodes that have values that are greatly 


ШОО "που апа опе Of those is directly under the load. 


Cee EOUSSINESQ PROBLEM 
The Boussinesq problem consists of one concentrated 
load normal to the surface of a semi-infinite solid as 


shown in Figure 23. 
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I Classical" Solution 
Inezelassieanlzsolmtion to the problem of a 
concentrated load normal to the surface of a semi-infinite 
solid was solved by J. Boussinesq [4] and yields the results 
shown in equations 20. The К СОТО Оеп: үре латте есг оп 


15 given in equation 21. 





ЈЕ = = (1-29) 25 - Eta V ?1-3r?z (92 422)797?] 
T T° 
E Εν, .’,.6Θ 6--575 
БЕРЕ τρ σα τσ) 
(20) 
σα = зу (1-2»){——% + rtz + z(e?) 
7 Т? г 
Ze - > р22(р2+52)72/2 
м = — [(14v)z?(r^422)793/2 4 2(1-52) (рожа? )71/23 (21) 


equations 20 and 21, E is Young's Modulus, and v is 
Poisson's ratio. The stresses are in cylindrical coordi- 
nates with the axis of symmetry being the line of action 
or che load. 

ο πι ος Formulation 

The TRISOP formulation took advantage of the double 

ΠΤ π᾿ πε problem, ana is shown in Figure оү» 
figure also shows the dimensions of the 3x3x3 mesh. The 
2x2x2, and 4x4x4 meshes were similarly constructed, and 
Figure 25 shows one face of these meshes to indicate the 


element dimensions. 
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Figure 24.  Boussinesq Problem TRISOP Formulation 
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Figure 25. Boussinesq Meshes 
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3. Results 

τε πο re O Dbypequatlon 2]. the deflection at the 
load is not bounded. After examining deflections at other 
points in the three meshes, it could not be determined if 
extrapolation to a converged solution would be valid Гог 
the meshes used. Figure 26 shows the solution of the 
Altec tions 1η the Z direction from r = 0 to r = 100, at 
EE oneescoft >, 10, and 20 inches from the surface using 
nodal point deflections from the 4x4x4 mesh. The TRISOP 
deflections follow the general contour of the classical 
solution, but are not as large. 

Figures 27, 28, 29, and 30 show contour graphs of 
sigma r, sigma 06, sigmaz and tau rz using classical and TRISOP 
data. A block twenty inches on a side is used for these 
graphs, because the area close to the load is of prime 
n σου, The data used in generating the TRISOP graphs 
fee strom interior nodal points in the mesh. Since TRISOP 
πο τος results in rectangular coordinates, and the 
For 2 mesoa,sommrıon is in cylindrical coordinates, it was 
necessary to transform stress data from the interior nodes 
КОШО осе се coordinates for meaningful comparison with 
μα ο ος σα solution. 

H. Conclusions 

The discrepancy between the classical solution for 
displacements and the TRISOP solution is caused by the 
ШОШ ОШ ο Ισ τος used. Іп спе real problem the deflec- 


tions at Z = 100 inches are not zero. If the boundary 
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Figure 26.  Boussinesq Deflections 
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Figure 29. Boussinesq Problem 
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пошао а, the base of the mesh could have modeled the 
expected deflections, it is believed that the deflections 
elsewhere in the mesh would have been nearer to the expected 
values. The contour graphs indicate that TRISOP generates 


Пошта in close agreement with the classical solution. 


0. PINCHED CYLINDER 

ИТ а о ес у оде consists of a thin shell cylinder 
with a concentrated radial load as shown in Figure 31. 
D. E. Hanson [2] solved this problem using TRISOP with four 
ο ος ροἱπς integration, and compared his results to a 
study mady be G. Cantin [7] using thin shell elements. The 
objective of this study was to determine if the three 
mensional solution would give comparable results. Hanson 
made runs with meshes up to a 1х10х10, and obtained 
discouraging results. After discovering that two Gauss 
о integration might yield better results, Hanson's 
1x2x2 and 1x4x4 meshes were rerun using two point integra- 
tion. Table 1 shows a tabulation of Hanson's and Cantin's 
Γ᾽ ἹπσιπατΠπςρ the reruns using two point integration. 
Figure 32 is a convergence study of the two meshes run with 
two point integration. The extrapolated result agrees with 


the fully converged solution given by Cantin. 
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ΡΤ 3l: 


Pinched Cylinder 
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10.35 inches 
4.906 inches 
0.09% inches 
100 pounds 





TABLE I 


Displacement of a Pinched Cylinder at the Applied Load 


Cantin Hanson 

NO. No. 

of Mesh Disp lacs, of Mesh Disp laca 

Eq. In.: Ба. IE 
36 Je -0.0921 193 хе хе -0.00177 
48 155, -0.1072 465 1x4x4 -0.00375 
60 1x4 -0.1099 045 1х6х6 -0.03002 
54 2х2 2020931 1593 1х8х8 ο 00 07 
72 2x4 -0.1085 2409 1x10x10 -0.10057 
90 2x4 -0.1113 

150 Hx! -0.1126 

204 6x6 -0.1137 

486 8x8 -0.1139 


726 021.0 -0.1139 


jH or σα τσαρτ πιο τοπ: Integration 


Νο. 

of Mesh Displae. 
Eq. μας 
153 202, -0.1043 
465 1x4x4 -0.1127 
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У. CONCLUSIONS AND RECOMMENDATIONS 


The analysis of the problems considered indicates that 
Or can be expected to give accurate results. However 
it was found in the simply supported beam and Boussinesq 
problems that data obtained for external surfaces were 
Шшамостгабе. Іс would therefore be wise to о πο 
data from surface nodes. If information near the surface 
of a solid is desired, a possible solution would be to have 
IS EA lanetro! elements at the surface, and use the 
data generated by the elements just below the surface. 
Another possible solution would be to modify the program 
MON compute stress and strain at points other than nodal 
τν: such as the Gauss integration points to see if 
Bester results are obtained. However, this would increase 
the complexity of the program, because the coordinates of 
шпенооїиоє to be used would have to be specified or calcu- 
lated and exhibited together with the stress values. 

ΙΙΙ; research LCDR E. Leonidas was 
making improvements to TRISOP [8] that included reducing 
the program run time, and generalizing the boundary condi- 
tions to make it possible to specify boundary displacements. 
This last improvement will make it possible to model 
problems such as the Boussinesq problem more accurately. 

πε Εμ ος ος ο ability οὗ receiving input mesh data 
ο πο τοι ο «πο, but this data is converted, and 


the output is given in rectangular coordinates. For problems 


29 





Имејл Boussinesqsproblem, output data in cylindrical 
coordinates would be a great help. TRISOP could be modified 
bo accomplish this, or possibly another similar program 
Шире Produced tor solve problems in cylindrical and 
spherical coordinates. 

The change from four point to two point integration 
τ Ξ; ο σσπαΙσεςρ time needed to solve a problem. 
This change also seems to give a more accurate solution as 
shown by the results of the pinched cylinder analysis. The 
Меше оно ше певавеа to solve large problems, however, can 
be excessive and is a very real limitation. Table ТТ shows 
the CPU time for various problems run with both CUB2 and 
CUB4 on an IBM 360 computer. 

neta rsest shortcoming of TRISOP is that the excessive 
amount of data produced is very time consuming and tedious 
to analyze. There are two possible solutions to this 
problem. The first would be to modify the program so that 
PeemouupuLy Would be placed on contour graphs Кор ен о а о ен 
Шшевергрорлеп маса спала solution is that if the ee deznor 
ОО ШУ or graphs other than those asked for are needed after 
reee program would have to be rerun. In the case 
“ПК Кре proplem, the cost in computer time could be 
excessive. What appears to be a far better solution would 
be to have the output placed on magnetic tape. The data 
could then be thoroughly analyzed using the computer with 


a minimum expenditure of computer time. 
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TABLE II 


ο FOR TEST PROBLEMS RUN 


Overall 
FSTF Merge Sole Time 
Simply Supported Beam (sec) (sec) (sec) Amin, sec) 
iix2 xc  CUB4 БЕЗ 120.77 Soon ie, 53258 
J> xO CUBI ΙΤ mms 539.26 20 06.58 
ШОЛ ХО ООВ! ОШО 353. 20 1116.13 JOE 
2x3x8 CUB! ЕШ 22 507,202 1357.63 Д7 44,98 
о сова 48.49 118,35 386.38 MEE SOUS 
Шх xo CUB2 ШЕТІ 58 534.60 [Sd 
О ООВ а D A БҶ DONE OE? 
О CUB БҮЙІ 59-20 1357.07 35 01.59 
Pinchea Disk 
1x4x4 CUBY ЭЛ ον 75.62 απο 177 ο δο 
ο CUB! 170.96 285.88 835.18 ΤΙ 6 
1хбхб СОВА 675.46 566.05 1694.84 H9 25.42 
1x8x8 CUBY DUE EM 38 4217.03 115 33.68 
1x4x4 CUB2 зе 173.04 502.34 Ш е с 
ІР ους» (fee 258.56 832.60 PO o 
I2 o CUB? πο 05656. 327 T6989.59 ΘΙ 
1x8x8 CUB2 КОШО 1133,12 4216.28 98 34.08 
Boussinesq Problem 
S2x x2 CUB2 25.19 90.13 249.43 G η 
Ол CUB2 ur E o == 1062.09 а 
4x4x4 CUB2 Ii 0655 60 4309.12 107 07.95 
Pinehea Cylinder 
1x4x4 CUBY ЕЕЕ 0 646.10 18 53.85 
IES UE? TS 497.40 12 42.00 
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